IGIS Tech Notes describe workflows and techniques for using geospatial science and technologies in research and extension. They are works in progress, and we welcome feedback and comments.

Background

The US Census Bureau reports most of their data by Census Tracts and Blocks, but they also maintain GIS data for many administrative boundaries which they use for other kinds of data summaries. You can download these geographies as GIS data.

In this Tech Note, we will download the boundaries of incorporated cities and towns in California, using the Census API. We will do this in R using the tigris package by Kyle Walker.


1. Load packages:

library(tigris)
library(dplyr)
library(sf)
library(leaflet)


Set Tigris cache option:

options(tigris_use_cache = TRUE)


2/. Get the Places Dataset

Incorporated cities and towns are part of the Places dataset, which we can download via places().

ca_places_sf <- places(state = "CA", cb = FALSE, progress_bar = FALSE) |> st_transform(4326)
## Retrieving data for the year 2021
dim(ca_places_sf)
## [1] 1611   17
head(ca_places_sf)
## Simple feature collection with 6 features and 16 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -120.6741 ymin: 33.62517 xmax: -117.6746 ymax: 40.45356
## Geodetic CRS:  WGS 84
##   STATEFP PLACEFP  PLACENS   GEOID         NAME          NAMELSAD LSAD CLASSFP
## 1      06   77364 02412017 0677364   Susanville   Susanville city   25      C1
## 2      06   02000 02409704 0602000      Anaheim      Anaheim city   25      C1
## 3      06   16532 02410239 0616532   Costa Mesa   Costa Mesa city   25      C1
## 4      06   17750 02410282 0617750      Cypress      Cypress city   25      C1
## 5      06   28000 02410556 0628000    Fullerton    Fullerton city   25      C1
## 6      06   29000 02410568 0629000 Garden Grove Garden Grove city   25      C1
##   PCICBSA PCINECTA MTFCC FUNCSTAT     ALAND  AWATER    INTPTLAT     INTPTLON
## 1       Y        N G4110        A  20507839  218688 +40.4336740 -120.6283062
## 2       Y        N G4110        A 130206498 1567025 +33.8555018 -117.7586572
## 3       Y        N G4110        A  40937896   24554 +33.6659055 -117.9123358
## 4       N        N G4110        A  17128118   22347 +33.8184768 -118.0383074
## 5       N        N G4110        A  58065841   32825 +33.8857156 -117.9280247
## 6       N        N G4110        A  46520862   31066 +33.7787816 -117.9604719
##                         geometry
## 1 MULTIPOLYGON (((-120.5261 4...
## 2 MULTIPOLYGON (((-118.0175 3...
## 3 MULTIPOLYGON (((-117.9546 3...
## 4 MULTIPOLYGON (((-118.0633 3...
## 5 MULTIPOLYGON (((-117.9854 3...
## 6 MULTIPOLYGON (((-118.0429 3...


Reading the docs, we can use the Legal/Statistical Area Description (LSAD) column to pull out the incorporated cities and towns. Let see what it contains:

ca_places_sf$LSAD |> table()
## 
##   25   43   57 
##  461   21 1129


Looking up the codes, we find:

25. City (suffix). Consolidated City, County or Equivalent Feature, County Subdivision, Economic Census Place, Incorporated Place

43. Town (suffix). County Subdivision, Economic Census Place, Incorporated Place

57. CDP (suffix). Census Designated Place, Economic Census Place (not incorporated)


Therefore, we want the Places where the LSAD is 25 (cities) or 43 (towns):

ca_cities_towns_sf <- ca_places_sf |> 
  filter(LSAD %in% c(25, 43))

dim(ca_cities_towns_sf)
## [1] 482  17


3. Plot the results

To plot the results, first download the state boundary:

cabnd_sf <- states(cb = TRUE) |> filter(GEOID == "06") |> st_transform(4326)
## Retrieving data for the year 2021


Make a version of the cities & towns layer that includes a column of color values:

random_cols <- c("#89C5DA", "#DA5724", "#74D944", "#CE50CA", "#3F4921", "#C0717C", "#CBD588", "#5F7FC7", "#673770", "#D3D93E", "#38333E", "#508578", "#D7C1B1", "#689030", "#AD6F3B", "#CD9BCD", "#D14285", "#6DDE88", "#652926", "#7FDCC0", "#C84248", "#8569D5", "#5E738F", "#D1A33D", "#8A7C64", "#599861")

ca_cities_towns_4map_sf <- ca_cities_towns_sf |> 
  mutate(col = sample(random_cols, nrow(ca_cities_towns_sf), replace = TRUE)) |> 
  select(NAMELSAD, col)


Plot them in leaflet:

leaflet(cabnd_sf) |> 
  addProviderTiles("CartoDB.Positron") |> 
  addPolygons(fill = FALSE, color = "black", weight = 5) |> 
  addPolygons(data = ca_cities_towns_4map_sf, color = "#777", weight = 1,
              fill = TRUE, fillColor = ~col, popup = ~NAMELSAD)


Looks good!

4. Export as CSV and GeoJSON

ca_cities_towns_sf |> 
  select(GEOID, NAME, LSAD, LAND_AREA_M2 = ALAND) |> 
  st_drop_geometry() |> 
  write.csv(file = "./outputs/ca_cities_towns_2021.csv", row.names = FALSE)


Save as GeoJSON:

ca_cities_towns_sf |> 
  select(GEOID, NAME, LSAD, LAND_AREA_M2 = ALAND) |> 
  st_write(dsn = "./outputs/ca_cities_towns_2021.geojson", delete_dsn = TRUE)
## Deleting source `./outputs/ca_cities_towns_2021.geojson' using driver `GeoJSON'
## Writing layer `ca_cities_towns_2021' to data source 
##   `./outputs/ca_cities_towns_2021.geojson' using driver `GeoJSON'
## Writing 482 features with 4 fields and geometry type Multi Polygon.


Done!

This work was supported by the USDA - National Institute of Food and Agriculture (Hatch Project 1015742; Powers).